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ABSTRACT 

Many  meteorological  forecast  applications  require  the  use  of  grids 
that  have  a  high  resolution  in  a  particular  area  of  interest,  while 
allowing  coarser  resolution  elsewhere.   Conventional  finite  difference 
models  often  use  nested  grids  to  this  end.   In  recent  years,  finite 
element  models  have  been  offered  as  an  alternative.   In  this  study,  the 
two-dimensional  advection  equation  with  diffusion  is  defined  over  a 
rectangular  domain.   The  Galerkin  technique  is  applied  to  linear  basis 
functions  on  triangular  elements.   The  model  is  tested  to  determine  the 
sensitivity  of  the  forecast  to  various  nodal  geometries.   Both  equi- 
lateral and  right  triangular  elements  are  tested.   It  is  found  that  the 
equilateral  arrangement  consistently  yields  a  superior  forecast.   Other 
tests  are  conducted  in  which  the  resolution  is  varied  smoothly  versus 
abruptly  over  the  domain.   The  smoothly  varying  case  gives  results  that 
are  dramatically  improved  over  the  abruptly  varying  case.   Among  the 
conclusions  is  the  fact  that,  for  a  given  maximum  resolution,  the  more 
slowly  and  smoothly  the  element  size  is  changed,  the  better  the  forecast 
obtained. 
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I.   INTRODUCTION 

Sub-Synoptic  scale  meteorological  elements  are  usually  not  repre- 
sented or  forecast  well  by  coarse  mesh  numerical  models.   This  situation 
may  be  improved  by  increasing  the  resolution  of  the  grid.   However,  a 
uniform  reduction  of  the  grid  mesh  requires  a  significant  increase  in 
computational  time  and  computer  storage.   Hence  this  is  a  limitation  on 
the  practical  resolution,  and  therefore  accuracy,  for  a  uniform  grid. 
The  conventional  solution  is  to  nest  grids  so  that  the  geographical  area 
or  meteorological  feature  of  interest  is  covered  by  a  fine  mesh,  and 
less  interesting  areas  or  features  are  left  with  a  coarse  mesh.   The 
advantage  of  this  technique  is  to  increase  accuracy  in  the  desired  area 
while  keeping  the  computational  storage  and  time  requirements  within 
reasonable  limits.   There  are,  however,  two  disadvantages  to  nested 
grids.   First,  the  abrupt  change  in  grid  size  at  the  boundaries  of  the 
fine  mesh  generates  noise  in  the  solution.   Second,  boundary  conditions 
for  the  fine  mesh  must  be  interpolated  from  the  coarse  mesh. 

A  better  technique  would  allow  the  grid  to  vary  smoothly  and  con- 
tinuously over  the  domain  rather  than  in  abrupt  jumps.   However,  such  a 
highly  variable  grid  greatly  complicates  a  finite  difference  numerical 
model.   As  an  alternative  to  such  a  finite  difference  technique,  consider 
the  Finite  Element  Method  (FEM) .   This  method  has  long  been  used  in 
mechanical  engineering  but  has  been  adapted  to  meteorology  only  during 
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the  last  decade.   Pioneering  work  in  meteorological  applications  was 
done  by  Cullen  (1973).   Staniforth  and  Mitchell  (1978)  and  Staniforth 
and  Daley  (1977)  have  devised  more  recent  forecast  models.   The  most 
recent  FEM  meteorological  model  at  the  Naval  Postgraduate  School  was 
written  by  Kelly  and  Williams  (1976)  and  is  the  precursor  to  this  study. 

The  great  attractions  of  the  FEM  are  the  more  accurate  phase  speeds 
[Haltiner  and  Williams  (1980)],  and  its  suitability  to  non-uniform  grids. 
The  method  involves  the  division  of  the  domain  into  a  number  of  elements 
and  then  defining  an  equation  for  each  node  of  each  element.   Thus, 
solution  requires  solving  a  large  number  of  linear  equations  simul- 
taneously.  However,  there  is  no  limitation  on  the  size  of  the  elements 
or  the  shape  of  the  domain.   The  method  is  easily  adapted  to  concentra- 
tion of  small  elements  in  areas  of  interest  and  assignment  of  larger 
elements  to  other  areas. 

The  purpose  of  this  investigation  is  to  determine  the  sensitivity  of 
a  FEM  model  to  different  grid  geometries.   The  model  used  is  a  two- 
dimensional  advection  model  with  diffusion.   It  was  deliberately  chosen 
because  of  its  simplicity.   By  keeping  the  model  simple,  all  variations 
in  the  results  can  be  attributed  to  differences  in  geometric  and  initial 
conditions  rather  than  inaccuracies  in  the  equation  formulation.   In 
addition  to  sensitivity  studies,  the  model  can  be  adapted  as  an  opera- 
tional predictor  of  any  scalar  field  such  as  cloudiness,  vorticity, 
temperature  or  precipitation. 
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The  Galerkin  technique  is  employed  to  transform  the  differential 
equation  for  the  dependent  variable  into  a  set  of  algebraic  equations 
Triangular  elements  are  used  with  the  area  coordinate  system.   Later 
chapters  will  show  that  these  techniques  greatly  simplify  the 
computations . 


13 


II.   THE  METHOD 

The  Galerkin  technique  will  be  explained  by  first  an  example,  after 
Crandall  (1956) ,  and  then  the  application  to  the  FEM,  after  Kelly  and 
Williams  (1976) . 

A.   EXAMPLE 

Consider  a  simple  differential  equation,  that  which  governs  the  cool- 
ing of  an  object 

x  =  1  t  =  0 

-  =  -x  t  >  0 

with  the  solution  defined  on  the  interval  0  <  t  <  1. 

The  critical  step  is  to  select  a  trial  family  of  approximate  solutions. 
(The  members  of  a  trial  family  are  often  called  basis  functions.)   For 
simplicity,  consider  the  second  order  trial-solutions 

2 
x  =  1  +  c  t  +  c  t  II-2 

where  x  is  the  approximate  solution. 

The  object  is  now  to  determine  the  coefficients  c,  and  c  .   Next,  form 

12 

the  residual,  R(t)  from  II-l 


dx 
R  =  —  +  x  II-3 

at 
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substituting  II-2  into  II-3 

R  =  1  +  c  (1  +  t)  +  c2  (2t  +  t2)  11-4 

Now  we  want  to  adjust  c  and  c   so  that  II-4  stays  close  to  zero  on  the 
interval  0  <  t  <  1. 

The  above  discussion  applies  to  any  weighted  residual  technique; 
however,  Galerkin's  technique  has  two  requirements: 

1.  The  weighted  averages   of  the  residual  must  vanish  over  the  interval. 

2.  The  weighting  functions   must  be  the  same  functions  of  t  as  were  used 

to  construct  the  trial  family  (in  II-2) .   In  this  case  those 

2 
weighting  functions  are  t  and  t  . 

That  is,  the  weighted  averages  are  formed  by  multiplying  the  residual 

by  the  weighting  functions,  integrating  the  product  over  the  interval 

and  setting  the  result  equal  to  zero. 


/ 


1 

t  R  dt  =  0 


o 
and  II-5 


1 


/ 


2 

t  R  dt  =  0 
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Substituting  II-4  into  II-5 


1  1 


/     t   R  dt   =     /      [t  +   c1 


(t+t2)    +   c2(2t2+t3)]dt 


=  —  +  —  c      +  —  c      =    0  II-6a 

2        6      1        12      2 


and 


/    t2R  dt  =     J      [t2+  C;L(t2+t3)  + 


c    C2t3+t4)]dt 


■  i +  n  ci +  To  c2  =  °  XI-6b 


Solving  II-6a  and  b  for  c  and  c  we  find: 


c  =  -.9143,  c2  =  0.2857 


So  that  the  approximate  solution  is 

x  =  1  -  0.9143  t  +  0.2857  t2 
The  exact  solution  to  Equation  II-l  is 


-t 
x  =  e 
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Figure  1  shows  the  plots  of  the  exact  solution  and  the  approximate  solu- 
tion.  Note  that  this  technique  yields  the  best  fit  only  on  the  desired 
interval  0  <  t  <  1  and  that  the  solutions  may  diverge  from  the  interval. 

In  general,  Galerkin's  technique  will  reduce  a  partial  differential 
equation  to  a  family  of  N  algebraic  equations,  where  N  is  the  number  of 

basis  functions.   In  this  example  Galerkin's  technique  has  reduced  II-l 

2 
to  two  equations,  I I -6a  and  II-6b,  in  the  two  basis  functions  t  and  t  . 

B.   APPLICATION  TO  THE  FEM 

The  Finite  Element  Method  divides  the  domain  into  discrete  segments 
called  elements  with  points  (called  nodes)  arranged  about  the  perimeter  of 
the  elements.   The  basis  functions  are  then  defined  locally.   These  basis 
functions  are  usually  low  order  polynomials  that  must  be  piece-wise 
continuous.   A  one-dimensional  example  is  Figure  2  wherein  the  domain 
(x-axis)  is  divided  into  four  elements  (line  segments)  W  through  Z.   All 
of  the  basis  functions  (A-C)  are  linear  rising  to  a  value  of  unity  over 
each  node  (points  2  through  4)  and  are  zero  elsewhere.   Now,  if  a  variable 
S  is  defined  over  node  3,  for  example,  then  it  may  be  approximated  by 


S  =  Y  V  +  v  V  +  Y„V  II-7 

'2  A    '3  B    '4  C 


where  y  ,  y  and  y  are  the  values  of  S  at  nodes  2,  3  and  4  respectively 
and  V  ,  V  ,  V  are  the  basis  functions  A,  B,  and  C  respectively. 

A     £3     \* 

Equation  II-7  may  be  substituted  into  any  equation  in  which  S  appears. 
An  equation  using  II-7  can  be  written  for  each  node,  by  multiplying  the 
basic  equation  by  the  basis  function  and  integrating  over  the  domain. 
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3.P 


Figure  1.   Solutions  to  cooling  object  example.   Approximate  solution 
is  the  curve  connecting  the  'x'  signs.   Exact  solution  is 
the  curve  connecting  the  '+'  signs. 
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/ 


Figure  2.   One-dimensional  example  of  linear  basis  functions, 


Figure  3.   A  domain  with  eight  nodes  and  eight  triangular  elements. 
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The  end  result  is  a  system  of  N  equations  in  N  unknowns,  where  N  is  the 
number  of  nodes.   Fortunately,  in  matrix  form,  these  equations  are  often 
tridiagonal  and  may  be  solved  with  great  economy. 

In  two  dimensions  the  basis  functions  are  somewhat  different,  but  the 
mathematics  is  identical.   In  Figure  3,  the  domain  has  been  divided  into 
eight  triangular  elements  with  eight  nodes.   Figure  4  shows  the  basis 
function  (outlined  in  heavy  black)  at  node  7  that  must  span  all  the 
elements  about  node  7.   Note  that  the  function  has  value  of  unity  at 
node  7  and  decreases  linearly  to  zero  at  all  surrounding  nodes.   Figure  5 
shows  another  basis  function,  but  this  one  at  node  5.   Each  node  has  a 
similar  basis  function  about  it. 

The  value  of  a  variable  S  may  be  approximated  at  any  node  i  by 


■£  yj  vi 


S  =  \    Y.   V.  II-3 

J 


where  j  ranges  over  all  the  nodes  connected  to  node  i  including  i  itself, 
Y.  is  the  value  of  S  at  node  j  and  V.  is  the  value,  at  the  j-th  node, 

d  : 

of  the  basis  function  of  the  i-th  node. 

As  in  the  one-dimensional  case,  an  equation  like  II-8  is  written 

for  each  node,  i,  then  multiplied  by  the  basis  function  and  integrated 

over  the  domain.   The  resultant  equation  set,  in  matrix  form,  contains 

an  n  by  n  square  matrix  where  n  is  the  number  of  nodes.   Although  this 

matrix  is  usually  not  tridiagonal,  careful  selection  of  a  nodal  numbering 

system  would  yield  a  matrix  that  is  strongly  banded  and  fairly  easily 

solvable.   Alternatively,  with  these  sparse  matricies,  iterative  methods 

converge  rapidly. 
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Figure  4.   The  linear  basis  function  associated  with  node  seven, 


Figure  5.   Another  linear  basis  function  on  the  same  domain. 
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III.   EQUATION  FORMULATION 

The  equation  of  interest  to  the  investigation  is  the  two-dimensional 
advection  equation  with  second  order  diffusion: 


3s     3s     3s      „2 

^-  +  u  -^-  +  v  -* K  V  S  =  0  III-l 

3t      ax  ay  h 


Define  an  inner  product  of  two  functions  f (x,y)  and  g(x,y)  as 

/   /  fg  dx  dy  5  <f,g> 
y  x 

Following  the  Galerkin  technique,  form  the  residue  and  find  the  weighted 
averages,  using  the  basis  functions  as  weighting  functions.  Neglect  the 
diffusion  term  at  this  stage. 


// 


y  x 


||  +  u  P.  +  v  P-]  V  dx  dy  =  0 

at  ox  ay 


or 


3s  3s  3s 

<-rr  ,  V>  +  <u  -s—  ,  V>  +  <v  -r—  ,  V>  =  0  III-2 

at  ax  ay 


where  V  is  the  weighting  function.   Although  the  domain  of  integration 
is  over  the  entire  range  of  x  and  y,  it  will  soon  be  shown  that  V  is  zero 
everywhere  except  locally  around  a  specific  node.   So  that  for  any  node, 
HI-2  is  non-zero  only  when  within  one  grid  increment  of  that  node. 
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The  variables  will  have  the  form 


S  =  Y.  V. 

3       3 


u  =  a.  V.  III-3 

J      J 


v  =  S.  V. 
3      3 


where  the  repeated  subscripts  indicate  summation  over  the  range  of  the 
subscript.   The  coefficient  is  a  function  of  time  only  and  the  basis 
function  is  a  function  of  space  only.   That  is 

S  =  y.  v.  =  Z     Y- (t)  V. (x,y) 
3      3        j    3  3 

etc. 

The  critical  step  is  the  requirement  that  the  weighting  function  in  III-2 
be  the  same  as  the  basis  function  in  III-3.   Substitution  of  III-3  into 
III-2  yields 

3v. 


2 

3j'i      kk'^dx   '   l 


<Y,v,,v,>  +  <akvkY  --,  v4> 


3v. 

+   <0,V  v.  T-*-  '    V.>  =  0  III-4 

k  k  ]  dy     i 


where  the  dot  indicates  a  time  derivative. 

Applying  first  the  Galerkin  technique  and  then  the  Gauss  Divergence 
Theorem  to  the  diffusion  term  yields 
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I  I  \  (V2S)V.  dxdy  =  \  /  /V  '(VS)V±   dxdy 
y  x  y  x 


=  k        I   /[V'(V.VS)  -  VS  •  VV.]  dxdy        III-5 

y  x 

=  IC  {    V.VS  •  n  dr  -   /  /  W.  •  VS  dxdy} 

y  x 

where  n  is  a  unit  vector  normal  to  the  domain  and  dr  is  the  differential 
distance  along  the  path  of  integration  on  the  perimeter  of  the  domain. 

as 

The  contour  integral  in  III-5  vanishes  because  -s—  cancels  out  due  to 

dx 

3S 
cyclic  continuity  and  ■»■—  =  0  because  there  is  no  flux  across  the  channel 

walls  at  the  north  and  south  edges  of  the  domain.   Using  III-3,  the  re- 
maining term  on  the  r.h.s.  of  III-5  can  be  written 

3v.       3v.      3V.       3V. 
-VW..VS>  -  -  V<si  ,  Y.  ^   ♦  <^  ,   Y.  ^1>]  in-6 

Combine  III-4  and  III-6  into  the  transformed  advection  equation  as 
follows 


3v.  3v. 

Y.<v.,v.>  +  y.{a  <v,  J-  ,  v.>  +  3,  <v,  ~— 2-  ,  v.> 


j  j   i     j   k  k  3x   '   i     k  k  3y   '   i 

3v.  3v 

i 

3y   '  3y 


3v.   3v.     3v.   3v. 
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Now  using  a  centered  difference  in  time  gives  the  forecast  equation 

3v.  3V. 

(Yn+1-Yn"1)<V.,V.>  =  -2At  Y.n(^n  <v    x-1  ,  v.>  +  Svn  <vv^-,v.> 

1  j     j     3   i         j        k    k  dx     i     k    k  tfy    i 


3v.        3v.  3v.        3v\ 

Kh[<^T  '  "3lT>  +  "17"  '  W 


+  K,  [<*r^   ,  ^-i>  +  <^  ,  r-^>]}  in-7 


when  n  is  the  time  step.   In  this  model  a  and  3   are  not  forecast 
quantities,  but  are  either  specified  for  each  time  step  or  are  constant. 
They  represent  the  prescribed  wind  field  which  is  advecting  the  scalar  S, 

This  equation  is  valid  for  each  node  i  and  as  such  is  a  matrix 
equation  of  the  form 

[A]  {x}  =  {b} 


where 


[A]  =  <  V.  ,  V.  >       with  dimensions   n  by  n 


{x}  =  (  Y •    -  Y •    )    with  dimensions   n  by  1 


{b}  =  The  right-hand  side  of  III-7 

with  dimensions   n  by  1 


and  n  is  the  number  of  nodes. 

Note  that  all  of  the  inner  products  are  independent  of  time  so  that 
they  need  only  be  calculated  once  and  placed  in  mass  storage  to  be 
accessed  at  each  time  step. 
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IV.   EVALUATION  OF  INNER  PRODUCTS 

Even  though  the  choice  of  elements  in  this  investigation  is  triangles 
with  nodes  at  the  vertices,  this  is  by  no  means  the  only  choice  possible. 
The  elements  could  be  triangles  with  six  nodes  each;  rectangles  with  six, 
eight  or  twelve  nodes;  or  any  other  polygon  or  curved-sided  element.   More 
complex  elements  would  allow  higher  orders  of  continuity.   Bathe  and 
Wilson  (.1976)  describe  the  technique  of  isoparametric  interpolation  as 
applied  to  the  FEM.   However,  these  more  complex  elements  require  even 
more  complex  formulations.   The  mathematics  can  get  out  of  hand  very 
quickly  requiring  many  more  terms  and  complicated  integrations.   Also 
Galerkin's  method  is  not  the  only  weighted  residual  technique  available. 
Crandall  (1956)  describes  the  collocation,  subdomain,  and  least  squares 
methods  as  well.   Lower  order  continuity  is  deemed  to  be  acceptable  in 
this  investigation  because  the  grid  resolution  is  fine  enough  so  that 
the  change  in  variable  is  assumed  to  be  linear  between  grid  points. 
This  assumption  allows  the  selection  of  triangular  elements  and  the 
Galerkin  method,  each  of  which  will  significantly  simplify  the  mathematics, 

Figure  6  shows  a  typical  element  with  a  point  p  described  in  area 

coordinates.   Lines  are  drawn  from  p  to  each  of  the  vertices  (x . ,  y., 

3         3 

j=l,2,3)  dividing  the  element  into  the  areas  A.  (j=l,2,3). 

Define:    L.  =  A. /A 
]    3 

where      A  =  total  element  area 
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3   x3,y- 


3»*3 


xi.yi 


P(x,y) 
or 
PCU.L2.L3) 


Figure   6.      A  typical  element  illustrating   the   area  coordinate 
system. 


AV 


-> 


Figure  7.   An  element  with  part  of  a  basis  function. 
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Now 


or,  in  matrix  form 


Ll  +  L2  +  L3  "  l 


x  =  Vl  +  L2X2  +  L3X3 


y  =  Vl  +  L2y2  +  L3y3 


L7! 


3  J 


or,  solving  for  the  L  vector 


1 
2A 


2A 


2A 


2A 


3  J 


IV- 1 


where 


ai  =  X3  "  X2 


bl  =  y2  "  y3 


a2  =  Xl  "  X3 


•  b2  =  y3  "  Yl 


a3  "  X2  "  Xl 


b3  =  Yl  "  y2 


So  that  any  one  of  the  equations  in  IV-1  is  really 


L:  - 1  +  A x  +  li y 


IV- 2 
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But  by  the  chain  rule 


^IL.    ..  tl 


3x   3x  3L.    2A  3l. 
J        ] 

and  IV- 3 

^     oL .  »      a .  ,>, 
d       j  9      j  d 

3y  "  "3y~  3l  .  "  2A  3L . 
J  D 


where  the  repeated  subscript  indicates  a  summation  from  one  to  three. 

Consider  Figure  7.   The  heavy  black  lines  outline  an  element 

divided  into  area  coordinates  defining  point  P,  and  only  A  is  labeled. 

The  hatched  triangle  (labeled  V  )  is  that  part  of  the  basis  function 

associated  with  node  2  that  lies  over  the  element.   There  are  two  more 

basis  functions  called  V  and  V  that  are  associated  with  nodes  1  and 

3  respectively.   Both  of  these  functions  also  have  linear  sections  over 

the  element.   However,  they  have  been  omitted  from  Figure  7  to  avoid 

clutter.   The  altitude  h   is  defined  as  unity.   The  altitude  h  is  the 

value  of  V  at  P.   Recall  L  =  A  /A.   As  P  moves  from  node  2  to  the 

opposite  side  of  the  element,  h  varies  linearly  from  1  to  0.   Following 

the  same  path  L  also  varies  linearly  from  1  to  0.   So  that  at  any 

point  P  on  the  element  L„  =  h  =  V"  .   In  general,  L.  =  V..   This  is  the 

2        2  j  j 

first  great  advantage  of  triangular  elements. 

So  now,  using  IV- 3 

3V.    b   dV.    b„  dV.    b„  dV. 
l    1   l   _2_  _jl    3   l 


3x    2A  3l    2A  3l    2A  3l 
12         3 


b,  3L.    b0  3L.    bn  3L. 

1    l     2    i     3    l 

2A  3L   +  2A  3L   +  2A  3L 
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but 


9L. 

SET =  1         xf       i  =  j 
=  0  if         i  ji  j 


So 


3v.       b.  3v.        a. 

^—^  =  r-         and   similarly  ■* —  =  rrr  IV-4 

dx  2A  dy     2A 


Note  that  these  derivatives  are  dependent  upon  the  individual  elements 
and  are  constant  with  time.   Therefore  they  need  only  be  calculated  once 
and  be  stored  as  part  of  the  inner  products. 

The  second  great  advantage  of  triangular  elements  and  linear  basis 
functions  is  the  general  area  integration  formula 


/ 


Lm  L, n  dA  = m!n!  rr;  2A  IV-5 

1   2       (m  +  n  +  2)1 


or,    rewritten  as  an  example  of  an   inner  product: 


2  2 '  1  •  A 

<L      ,    L  >  -  =-±-  2A  =  — 
1    '       2  5!  30 


Now  the  inner  products  for  equation  III-7  may  be  rewritten  with  the 
aid  of  IV-4  and  IV-5: 


30 


1 '  1 '  A 

<V.,    V.>   =    <L.,    L.>    -   -~ -  2A  =   — 
j         1  j         l  4!  12 


if   i    *    j 


2.  2    1      _.  A 

<L.    >  =  — r-r-  2A  =  — - 

l  4!  6 


if   i   =    j 


3V.  b. 

k    3x  i  2A  k'       l  2A   12         24 


b.     n         b. 


if  k/i 


3v.      9v. 

3x       3x 


b.    b. 
2A   2A 


3v.      3v. 


3y   '    3y 


JL>  = 


a.    a . 

2A   2A 


b. 
=      3 

12 


if  k  =   i 


<V. 


3v.  a. 

k   3y    '      i  24 


if  k  ^   i 


a. 

=  _i 

12     . 


if  k   =   1 
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V.   INITIAL  AND  BOUNDARY  CONDITIONS 

Simple  initial  conditions  are  used  for  the  scalar  field  S  and  the 
wind  components  u  and  v.   The  S  field  consists  of  a  half  sine  wave  in 
the  y  direction  and  multiple  waves  in  the  x  direction.   These  multiple 
waves  are  the  mechanism  used  to  generate  smaller  scale  features.  The 
—  term  in  the  u  equation  causes  the  area  of  fastest  flow  to  be  coinci- 
dent with  the  area  of  highest  spatial  resolution. 


,2nTTx.   .  4  ,TTv 
S  =  A  sxn  (— — )  sin  (— )  V-l 

L  W 

—  27Tv    TT 

u  =  U  +  B  cos  Or—  -  £)  V-2 


v  =  V  =  0 


where 

A  =  arbitrary  amplitude  =  400  meters 

n  =  wave  numbers        =  1 , . . . ,  6 

L  =  channel  length      =  2400  km 

W  =  channel  width       =  2400  km 

U  =  mean  zonal  wind     =  20  m/sec 

B  =  U  x  R;  R  is  the  ratio  of  the  perturbation 
wind  to  the  mean  zonal  wind 

Figures  8  and  9  depict  the  initial  fields  of  S  and  u  respectively. 
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x   km 


Figure  8.   Initial  field  for  the  scalar  S. 
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Figure  9.   Initial  u  field  with,  wind  ratio  R  =  0.4. 
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Boundary  conditions  are  also  handled  simply.   Cyclic  continuity  is 
assured  in  the  x-direction  by  selecting  the  proper  node  convectivity. 
Two  techniques  are  used  to  control  flow  across  the  channel  walls  at 
the  north  and  south  sides.   First,  the  exponent  on  the  y-component  of 
the  initial  wave  is  high  enough  to  keep  the  value  of  S  close  to  zero 
in  the  vicinity  of  the  walls.   Second,  the  solution  method  is  iterative 
so  that  the  values  of  the  nodes  one  row  in  from  the  channel  walls  are 
set  equal  to  the  nodes  on  the  walls  at  each  iteration.   The  result  is 
that  diffusion  from  the  edges  into  the  interior  is  disallowed. 
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VI.   COMPUTATIONAL  TECHNIOUES 

A.   NUMBERING  SCHEMES 

As  described  in  Chapter  IV,  the  inner  products  may  be  easily 
evaluated  on  the  individual  element  level.   To  facilitate  this  evalua- 
tion, a  local  numbering  system  is  required  for  each  element.   An  array 
called  IEL,  dimensioned  N  by  3,  contains  the  identification  number  of 
the  nodes  on  the  three  verticies  for  each  of  N  elements.   For  each 
element  these  node  numbers  must  be  stored  sequentially  in  a  positive 
sense,  that  is  counterclockwise.   The  node  with  which  to  start  number- 
ing, however,  is  arbitrary. 

In  addition  to  the  local  numbering  system,  a  "mass"  matrix  is 
needed.   The  mass  matrix  is  a  matrix  of  coefficients  whose  rows  are 
the  equations  of  the  system  to  be  solved.   There  is  one  equation  (row) 
for  each  node.   Each  equation  will  have  a  term  (column)  for  each  node. 
A  non-zero  entry  in  the  mass  matrix  indicates  connectivity.   That  is, 
if  entry  (I, J)  is  non-zero  then  node  I  is  connected  to  node  J.   Two 
nodes  are  connected  if  they  are  both  vertices  of  the  same  triangular 
element.   In  this  model,  each  node  is  connected  to  a  minimum  of  two  and 
a  maximum  of  six  or  eight  other  nodes  depending  upon  model  version. 
The  model  uses  an  array  called  MAME  dimensioned  N  by  7  (or  9)  that  con- 
tains the  numbers  of  all  nodes  connected  to  any  specific  node,  includ- 
ing itself,  for  each  of  N  nodes.   MAME  will  contain  from  3  to  7  (or  9) 
non-zero  entries  in  each  row. 
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To  save  storage  space,  MAME  is  next  compacted  into  a  one  dimensional 
vector  NAME  which  contains  only  the  non-zero  entires  of  MAME.   The 
utility  vectors  ISTART(N)  and  NUM(N)  are  also  constructed  in  order  to 
decode  NAME.   ISTART(N)  contains  the  index  in  NAME  of  the  starting 
position  of  the  continuity  of  any  node  N.   NUM(N)  contains  the  length 
of  the  connectivity  of  node  N  within  NAME.   As  an  example  of  how  to  re- 
trieve the  connectivity  of  a  node,  say  node  14,  from  NAME:  Let  J  = 
ISTART(14)  and  let  K  =  NUM(14).   Then  NAME (J)  through  NAME(J+K-1)  contain 
the  connectivity  for  node  14. 

Consider  the  matrix  equation  [A]  {x}  =  {b}  from  Chapter  III.   It 
should  be  clear  that  [A]  is  the  mass  matrix  mentioned  above.   It  is  a 
very  sparse  matrix  containing  the  inner  products  <V.,  V.>  for  row  j  and 
column  i.   In  this  model  [A]  is  compacted  into  the  vector  {h}  using  NAME 
as  a  lookup  table.   The  result  is  that  there  is  no  longer  a  real  matrix 
equation,  but  rather  just  the  product  of  vectors 

{H>  {x}  =  {b}  VI-1 

B.   SOLUTION  TECHNIQUE 

An  advantage  of  economy  of  space  has  been  gained,  but  VI-1  now  re- 
quires an  iterative  solution.   The  solution  chosen  is  a  standard  Gauss- 
Seidel  technique  which  requires  a  reversal  of  direction  of  iteration 
every  pass.   An  average  of  16  passes  are  needed  for  each  time  step.   The 
solution  is  considered  to  have  converged  to  its  final  value  when: 
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x.   -  x.  , 

X. 

1 


for  every  node  i. 

with:    I     =  iteration  number 

x.  =  solution  for  node  i  in  x 

1 

More  efficient  solution  schemes  are  available,  but  this  one  was 
selected  because  of  its  simplicity.   Minimizing  computation  times  is 
not  of  high  priority  to  this  investigation  for  two  reasons:  first,  the 
recently  installed  IBM  3033  is  a  fast  machine,  and  secondly,  since  an 
advection  model  doesn't  contain  gravity  waves,  stability  can  be  main- 
tained with  very  large  time  steps.   As  an  example,  a  twelve  hour 
forecast  with  a  grid  size  of  200  km  and  a  time  step  of  one  hour  requires 
less  than  eight  seconds  of  CPU  time.   Even  the  largest  version  of  the 
model  with  longer  forecast  length,  higher  resolution  and  linkage  to 
the  plotting  software  takes  less  than  eight  minutes. 

C.   ROBERT  FILTER 

Haltiner  and  Williams  (1980)  discuss  the  advantages  of  time  averag- 
ing.  Averaging  operators  act  as  low-pass  filters  that  tend  to  remove 
high  frequency  waves  while  having  little  effect  on  long-period  waves. 
Even  very  weak  filters  will  remove  high  frequency  noise  if  used  at  every 
time  step.   In  this  model,  the  choice  of  filters  is  that  of  Robert  (1966) 
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First,  assume  that  an  average  field  S    ,  already  exists  at  time 

level  n-1  as  well  as  the  unaveraged  field  S  at  time  level  n.   Then 

n 

the  model  uses  its  predictor  equation  to  determine  the  solution  vector 
{x}.   For  each  node,  an  unaveraged  predicted  value  is  computed  as 


S   ,  =  S   ,  +  X 

n+1     n-1 

Next,  the  corrected  values  at  time  n  are  determined  from 

?  =S   +Y(S   ,  -  2S   +  i"  , )  VI-2 

n    n       n+1     n    n-1 

where  y  ^s  the   averaging  coefficient.   Finally,  S   is  stored  in  place 

of  S   ,  and  the  model  continues  on  to  the  next  step. 
n-1 

The  value  of  y  must  be  carefully  chosen  as  it  does  effect  the  com- 
putational stability.   As  y  increases,  the  maximum  time  step  decreases. 
Haltiner  and  Williams  prefer  a  y  less  than  0.25  in  order  to  permit  a 
reasonable  time  step.   Additionally,  it  should  be  noted  that  a  large 
value  of  y  will  begin  to  dampen  waves  in  the  meteorological  frequency 
range.   On  the  other  hand,  Robert  used  a  value  of  y  =  0.01  applied  over 
a  total  forecast  length  of  many  days.   This  very  weak  filter  was  all 
that  was  needed  for  noise  suppression. 

For  the  model  in  this  investigation  y  is  an  input  parameter.   For 
one  set  of  tests  in  Chapter  X,  the  value  of  y  is  varied  from  0.01  to 
0.04. 
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D.   DETERMINATION  OF  TIME  STEP 

Cullen  (1973)  calculated  the  stability  criterion  for  the  two-dimen- 
sional problem  as: 


,     Ax 
At  <  VI-3 


:i/6 


the  minimum  grid  distance,  Ax,  varies  from  200  km  to  about  40  km.   The 
propagation  speed,  c,  varies  from  40  m/sec  to  0.   Accordingly,  the  time 
step  has  been  chosen  as  60  minutes  for  the  coarse  mesh  versions,  and 
ten  or  five  minutes  in  the  finer  mesh  versions. 
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VII.   NODAL  GEOMETRIES 

Once  the  decision  has  been  made  to  use  triangular  elements,  the  next 
choice  is  the  type  of  triangle.   This  investigation  concerns  both  right 
and  equilateral  triangles. 

A.   RIGHT  TRIANGLES 

Figure  10  is  an  example  of  the  right  triangle  arrangement  used  in 
this  model  for  a  uniform  grid.   The  domain  is  divided  into  a  series  of 
rectangles,  and  then  each  rectangle  is  bisected  by  a  diagonal  to  form  a 
pair  of  triangles.   Kelly  and  Williams  found  that  significant  bias  was 
introduced  into  the  model  if  all  of  the  diagonals  sloped  in  the  same 
direction.   This  model  overcomes  that  kind  of  bias  by  alternating  the 
slope  of  the  diagonals  from  rectangle  to  rectangle  both  horizontally  and 
vertically.   Consequently,  the  connectivity  of  the  nodes  will  vary  from 
four  to  nine  total  connections,  and  each  node  connects  with  itself. 

Two  methods  are  used  to  vary  the  grid  resolutions  with  right 
triangles.   In  the  first  method  a  FORTRAN  DATA  statement  specifies  co- 
efficients used  in  the  calculation  of  the  nodal  coordinates.   This  method 
is  very  simple  and  has  the  advantage  of  keeping  all  the  nodes  in  a 
rectangular  arrangement.   Rectangularly  arranged  nodes  lend  themselves 
well  to  harmonic  analysis  with  some  interpolation.   The  grid  may  be 
varied  in  either  one  or  both  the  x  and  y-directions.   Figure  11  is  an 
example  of  the  nodal  arrangement  by  this  method  with  variation  in  both 
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Figure  10.  The  domain  divided  into  right  triangular  elements.  Only 
the  elements  in  the  lower-left  corner  are  drawn,  but  the 
pattern  is  repetitive  over  the  entire  domain.  Each  node 
is  marked  with  an  'X'. 
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Figure  11.   Same  as  figure  10  with  non-uniform  elements. 
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directions.   The  disadvantage  here  is  that  the  resolution  must  be  kept 
relatively  fine  in  parts  of  the  domain  away  from  the  area  of  interest 
(the  center) .   Optimumly,  the  resolution  should  be  fine  about  the  area 
of  interest  and  more  coarse  elsewhere,  hence  the  next  method. 

In  the  second  method,  the  domain  is  divided  into  quadrants  with  the 
origin  at  the  center  of  the  domain.   The  nodes  around  the  outer  boundary 
of  the  first  quadrant  are  specified  as  having  the  same  coordinates  as 
the  uniform  grid.   But  the  nodes  along  the  abscissa  and  ordinate  are 
compressed  toward  the  origin  by  the  use  of  another  DATA  statement.   Con- 
sider four  nodes  each  with  coordinates  (x.,  y.),  i  =  1...4  .   Node  1  is 

1   1 

on  the  ordinate,  node  2  on  the  right  boundary,  node  3  on  the  abscissa, 
and  node  4  on  the  top  boundary.   The  line  from  node  1  to  2  is  defined 
by 


y  -  yx    y2  -  Yj 

= VII-1 

x   x1  x2    x1 


And  the  line  connecting  nodes  3  and  node  4  is 


y  "  y3      Y4  I  Y3 
x  "  x3      x4  "  x3 


VII-2 


Equations  VII-1  and  VII-2  are  then  solved  simultaneously  to  calculate  the 
coordinates  of  the  node  located  at  the  intersection  of  the  two  lines.   By 
selecting  the  proper  pairs  of  nodes  along  the  boundaries  and  axes,  the 
coordinates  of  all  the  nodes  on  the  interior  of  the  first  quadrant  are 
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found  by  VII-1  and  VII-2.   The  nodes  about  the  origin  are  then  slightly 
adjusted  to  allow  local  uniformity.   Finally  the  first  quadrant  is  re- 
flected about  the  axes  to  form  the  other  three  quadrants  in  the  domain. 
Figure  12  is  an  example  of  such  an  arrangement. 

The  advantage  of  this  method  is  that  the  resolution  is  fine  only  in 
the  area  of  interest,  and  changes  smoothly  in  all  directions  away  from 
that  area.   The  disadvantages  are  the  complexity  and  the  non-rectangu- 
larity  of  the  nodes. 

B.   EQUILATERAL  TRIANGLES 

Preliminary  results  of  the  model  indicated  the  presence  of  a  degree 
of  noise  when  using  right  triangles.   Cullen  (personal  conversation) 
suggested  the  use  of  equilateral  triangles  as  an  alternative  element 
arrangement.   Figure  13  is  an  example  of  such  an  arrangement.   The 
right  triangles  along  the  sides  of  the  domain  are  in  reality  halves  of 
equilateral  triangles.   Because  of  cyclic  continuity,  the  domain  actually 
"wraps  around"  so  that  all  of  the  elements  are  equilateral.   Note  that 
the  maximum  connectivity  is  now  only  seven  and  that  the  domain  is  no 
longer  square  but  rectangular. 

Variation  of  resolution  is  now  somewhat  more  complicated  and  is 
accomplished  by  a  transformation  of  coordinates.   The  desired  end  is  a 
smooth  and  easily  controllable  stretehage  and  shrinkage  of  the  coordi- 
nate axes.   In  the  x-direction,  the  transformation  is 

X  =  x  +  a  cos  kx  VII-3 
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Figure  12.   An  alternate  method  of  varying  resolution, 
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Figure  13.   Same  as  figure  10  except  with  equilateral  triangular 
elements. 
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where  X  =  transformed  coordinate 

x  =  original  coordinate 
a  =  coefficient  to  be  determined 

k  =  T 

L  =  channel  length 


3x 
the  map  factor  -*—  is  defined  by 
dx 


3X       u 

•*—  =  1  -  ak  sin  kx 

dx 


whose  maximum  and  minimum  values  are 


&  =  1  +  ak  ;  &    .   =  1  -  ak 

dx  max  dx  mxn 


The  ratio,  r  ,  of  maximum  stretch  to  minimum  shrink  is 

1  +  ak 
ri   1  -  ak 


or,  solving  for  a 


To  keep  the  transformed  coordinate  system  from  folding  back  on  itself, 
a  must  be  positive.   That  is,  r  ,  must  be  equal  to  or  greater  than  one. 
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Variation  of  resolution  is  accomplished  by  selecting  an  appropriate 
value  for  r  and  then  substituting  equation  VII-4  into  VII- 3.   Similar 
expressions  can  be  derived  for  the  transformation  of  the  y-coordinate 
by  the  use  of  the  ratio  r  . 

The  advantage  of  this  method  is  the  extremely  sensitive  control  upon 
resolution  afforded  by  the  ratios  r  and  r  .   Once  programmed,  this  method 
is  much  easier  to  use  than  the  more  cumbersome  DATA  statements  of  the 
previous  methods.   The  major  disadvantage  is  that  the  transformed  system 
is  no  longer  made  up  of  equilateral  triangles  if  one  of  the  ratios  is 
different  than  one.   Figure  14  is  the  nodal  arrangement  with  r  =  5  and 
r2=3. 
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Figure   14.      Same   as   figure 


13  but  r,  =  5,  r_  =  3, 
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Figure  15.   Initial  S  field  for  the  diffusion  tests, 
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VIII.   FORCING  TERM 

As  part  of  an  effort  to  improve  the  testing  of  a  numerical  scheme, 
it  is  often  beneficial  to  force  the  scheme  with  the  desired  solution, 
if  known.   Since  this  model  is  a  simple  advection  scheme,  the  output  will 
look  very  similar  to  the  initial  condition. 

Assume  the  solution  is  of  the  form 

S  =  cos  [u(X  -  At)  +  6] 


where 

27T 
y  =  Ix"' 

A  =  phase  speed, 

6  =  "  2  ' 

X  =  transformation  of  x  coordinate 


If  v  =  0,  then  the  forced  advection  equation  is 


■~  +  u  ~  =  F(x,t)  VIII-1 

dt      dx 


but 


3S  TT 

—  =  uA  sin[y(.X  -  At)  -  — ]  VIII-2 

dt  2 

and  • 

IS.  =  -u  |5  sin[y(X  -  At)  -  ~]  VIII-3 

dx      dx  2 
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Substituting  VIII-2   and  VIII-3   into  VIII-1 


u(A  -  u  ~)    sin[u(X  -   At)    -  |]    =  F(x,t) 


which   is   equivalent   to 


^  y  IT  TT 

y(A   -  u  -r— )  [sin(yx  -  -)    cos  (yAt)  -  cos  (yX  -  — )    sin(Ayt)]    =  F(x,t) 
dx  2  2 


applying  Galerkin's  technique  this  equation  eventually  reduces  to 


where 


(f.  -  g.)  <V.,  V.>  =  F(x,t)  VIII-4 

3    3    3   3. 


^X  TT 

f .  =  li(A  -  u  -r-)  sin(yX  -  — )  cos  yAt 
3  dx  2 


VIII-5 


g.  =  y(A  -  u  -r— )  cos  (yx  -  — )  sin  yAt 


Now,  VIII-4  should  be  included  on  the  r.h.s.  of  the  predictor  equation 
III-7.   Notice  that  both  of  the  equations  VIII-5  contain  a  time  depen- 
dent part.   The  technique  now  is  to  calculate  the  time  independent 
coefficients  of  VIII-5  at  the  onset  and  store  them  for  the  duration  of 
program  execution.   During  the  forecast  sequence,  at  each  time  step, 
the  time  dependent  parts  of  VIII-5  are  calculated,  multiplied  by  the 
applicable  pre-stored  coefficients  and  then  assembled  into  the  mass 
matrix  as  indicated  by  VIII-4. 
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This  forcing  term  will  not  be  used  in  most  of  the  model  runs  as  a 

standard  part  of  the  predictor  equation,  but  will  be  the  subject  of  a 

specific  set  of  tests  designed  to  determine  its  contribution  to 
forecast  accuracy. 
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IX.   ANALYTIC  SOLUTION 

Consider  the  one  dimensional  equivalent  of  equation  III-l  with  no 
diffusion: 

|{U  u(x)  |£-  0  IX-1 

3t        9x 

where  flow  is  zonal  only  and  the  velocity  may  vary  spatially.   Recall 
from  Chapter  V  that  the  zonal  wind  equation  is 

u(x)  =  U  +  B  cos(kx  -  -)  IX- 2 

where  k  =  — —  ,  L  =  zonal  wave  length. 

L 

In  a  simple  one  dimensional  advection  equation  for  any  time  t,  the 
following  holds: 

S(xft)  -  S(x  ,  0) 

where  x   is  some  point  upstream  in  the  initial  field.   Now  if  x   can  be 
determined  as  a  function  of  the  current  location  and  the  elapsed  time, 
then  the  analytic  solution  for  the  forecast  at  that  time  would  be 
achieved.   That  is, 


xQ  =  F(x,t) 
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dx 
Since  u  =  —  ,  then  from  IX-2 
at 


0    —  TT. 

— 7-  =  U  +  B  cos  (kx.  -  -) 
at  U    2. 


Integrating 


x  c 

I  ^ W 

7  U  +  B  cos(kxQ  -  j)    J 
X0 


kdt 


With  the  help  of  a  Table  of  Integrals,  this  expression  is  evaluated  as 


.           .  (U-B)tan  T(kx  -  ?-) 
(tan"1  [ r2— 2-] 


(U2-B2)1/2  (U2-B2)1/2 


,  (U-B)tan  |-(kxn  -  \) 


Let 


U  -  B 
r  = 


CU2  -  B2)1/2 


So  that 


tan"1  [r  tan  |(kx  -  |)  ] -tan"1  [r  tanj(kxQ-|)]  =  j(U2-bV/2  kt 
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Solving  for  x 


^    2     -1   rl        -1,      1„    TT,    1,-2    2  1/2    , , 
x„  =  —  +  r  tan    I—  tan  [tan   (rtan-(kx--)  -  -(U  -B  )     kt]  }  IX- 3 
0    2k   k         r  2      2     2 


The  analytic  solution  at  time  t  for  any  node  at  coordinates  (x,y) 
is  found  by  substituting  x  and  t  into  IX-3  to  obtain  x  .   Then  x  and 
y  are  substituted  into  equation  V-l.   This  process  requires  each  value 
of  x  to  be  operated  upon  by  five  trigonometric  functions,  each  of  which 
further  compounds  truncation  error.   As  a  result,  this  subroutine  must 
be  run  in  double  precision.   Even  so,  there  is  still  some  small  error 
in  the  analytic  solution  field.   This  investigation  is  concerned  with 
the  relative  errors  of  various  nodal  arrangements  and  is  not  intended 
for  comparison  studies  between  finite  element  and  other  methods.   Since 
the  errors  in  all  test  cases  are  relative  to  the  same  solution,  the 
small  error  in  the  solution  is  not  considered  significant. 


53 


X.   RESULTS 

The  model  contains  eight  variables:   element  type  (equilateral  and 
right  triangles) ,  east-west  resolution,  north-south  resolution,  wave 
number,  diffusion,  Robert  filter,  uniformity  of  flow,  and  number  of  nodes 
on  the  domain.   Additionally,  the  time  step  may  be  varied  and  the 
initial  conditions  may  be  changed  so  that  sharper  waves  are  created. 
If  all  of  these  variables  were  tested  throughout  their  meteorological 
range,  and  if  all  the  interactions  between  the  variables  were  investi- 
gated, the  number  of  computer  runs  required  would  be  well  over  ten 
thousand.   Since  that  many  runs  is  not  practical,  decisions  have  to  be 
made  concerning  the  scope  of  this  investigation.   There  seems  to  be  one 
of  two  ways  to  proceed: 

1.  One  or  two  questions  could  be  answered  rather  exactly  with  an 
attempt  to  quantify  the  relationship  between  only  a  couple  of  variables. 

2.  General  answers  can  be  attempted  for  many  more  questions  with 
the  aim  to  set  qualitative  guidelines  for  future  finite  element  studies. 

In  the  first  alternative  above,  quantitative  relationships  would  be 
of  little  value  unless  they  could  be  extended  to  FEM  models  as  a  class. 
Proof  that  such  an  extension  is  valid  is  clearly  beyond  this  investiga- 
tion.  As  a  result,  the  second  alternative  is  chosen  and  only  rather 
general  questions  are  qualitatively  answered. 

The  model  is  run  a  total  of  137  times.   The  runs  are  compared  and 
contrasted  in  a  number  of  test  cases.   Several  test  cases  comprise  the 
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investigation  of  each  question.   Two  types  of  statistics  are  calculated 
for  each  run.   First  of  all,  a  harmonic  analysis  is  accomplished  using 
nodes  along  each  latitude  circle,  for  the  initial  condition,  the  forecast 
field  and  the  analytic  solution  field.   Since  the  analysis  requires 
equally  spaced  values,  there  is  an  interpolator  that  calculates  values 
at  specific  points.   The  linearity  of  the  basis  functions  lends  itself 
very  well  to  linear  interpolation.   One  subroutine  determines  in  which 
element  any  particular  interpolation  point  lies.   Another  routine  uses 
the  coordinates  of  that  point  and  the  values  at  the  nodes  of  the  proper 
element  to  calculate  the  interpolated  value.   In  the  runs  that  include 
the  forcing  term  (see  Chapter  VTI) ,  the  interpolater  is  disconnected  and 
the  raw  nodal  values  are  used.   This  harmonic  analysis  generated  ampli- 
tude and  phase  information  for  all  possible  wave  numbers.   In  the  follow- 
ing paragraphs  the  term"phase  speed"  denotes  the  ratio  of  the  phase 
shift  of  the  forecast  wave  to  the  phase  shift  of  the  wave  in  the  analytic 
solution. 

The  second  type  of  statistics  generated  considers  the  forecast  and 
solution  values  as  an  ordered  pair  for  each  node.   A  root-mean-square- 
error  (RMSE) ,  correlation  and  bias  are  calculated  for  this  set  of  pairs. 
By-products  of  these  statistics  that  are  not  referred  to  hereafter  are 
the  applicable  means  and  standard  deviations. 

The  following  sections  are  organized  such  that  each  one  addresses 
an  individual  problem,  question  or  technique.   To  this  end,  within  each 
section,  one  or  more  of  the  eight  variables  listed  above  is  varied  from 
a  standard  or  nominal  model  configuration.   Changes  in  the  generated 
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statistics  are  then  attributed  to  the  influence  of  the  altered  variables. 
The  nominal  model  configuration  follows:   equilateral  elements,  uniform 
grid  (r  =  r  =  1,  see  Chapter  VII),  wave  number  one,  no  diffusion, 
Robert  filter,  y  =  0.1,  wind  perturbation  ratio  R  =  0.0  (see  Chapter  V), 
number  of  nodes  =  13  X  12.   The  standard  forecast  length  is  48  hours, 
during  that  period  any  feature,  moving  with  the  mean  flow,  should  move 
3456  km  or  around  the  domain  1.44  times.   Section  A  is  a  general  statement 
based  upon  nearly  all  of  the  runs.   It  includes  changes  made  to  all  of 
the  variables.   In  all  of  the  remaining  sections,  the  departure  from  the 
nominal  configuration  is  annotated. 

A.   ACCURACY  VS.  RESOLUTION 

For  this  test,  six  runs  are  discarded  because  they  are  deliberately 
constructed  to  test  unique  factors  discussed  later.   In  the  other  131 
runs,  the  forecast  values  stay  within  reasonable  limits  for  all  combina- 
tions of  resolution.   In  general,  the  accuracy  over  the  domain  as  a  whole 
decreases  as  the  stretch  of  the  grid  resolution  increases,  but  the 
statistics  all  remain  bounded.   The  correlations  are  all  above  0.96,  the 
phase  speeds  are  within  four  percent  and  the  RMSE  is  always  less  than 

eight  percent  of  the  range  of  the  amplitude.   Typical  values  are  cor- 
relation 0.99,  phase  speed  1.008  and  RMSE  one  percent  of  range. 

This  result  is  very  encouraging  in  that  no  cases  have  to  be  dis- 
regarded simply  because  they  can  not  be  explained.   The  cases  are 
consistent  with  one  another  and  exceptions  are  few. 
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B.   EFFECT  OF  DIFFUSION 

A  particularly  noisy  case  is  used  for  this  test.   The  domain  length 
scale  has  been  doubled  so  that  the  non-uniform  wind  is  more  obvious. 
With  resolution  of  r  =1/  r  =  2  wind  perturbation  ratio  R  =  0.4,  and 
13  X  24  nodes,  the  diffusion  values  are  incrementally  increased  from 
zero.   A  total  of  eight  runs  are  made.   Figure  15  is  the  initial  condi- 
tion for  this  test  and  Figure  16  is  the  48  hour  forecast  with  no  diffu- 
sion.  As  diffusion  is  increased  the  reduction  in  RMSE  is  dramatic  and 
the  correlation  increases,  but  the  wave  amplitude  diminishes  and  the 
maximum  time  step  must  be  decreased  to  insure  stability.   However,  a 
critical  point  is  eventually  reached  at  which  the  amplitude  is  reduced 
so  much  that  the  RMSE  starts  to  climb  again  and  correlation  falls. 
Figure  17  is  the  48  hour  forecast  with  the  value  of  diffusion  that 
nearly  minimizes  the  RMSE.   Note  the  absence  of  noise  but  the  obvious 
loss  of  amplitude.   Compare  this  with  Figure  18,  the  analytic  solution. 
Notice  that  the  contour  packing  is  downstream  of  the  high  in  Figures  16 
and  18,  but  upstream  of  the  high  in  Figure  17.   This  reversal  of  packing 
is  a  by-product  of  the  manner  in  which  the  model  applies  diffusion. 
Cullen  (personal  communication)  suggests  that  diffusion  be  increased  in 
the  areas  of  fastest  flow  and  decreased  where  the  flow  is  slower. 
Following  this  advice,  the  model  applies  diffusion  proportional  to  the 
windspeed  at  each  node.   The  effect  with  non-uniform  flow  is  that  diffu- 
sion is  also  non-uniform,  hence  the  slight  alteration  in  the  forecast 
field.   Diffusion  appears  to  be  a  very  effective  noise  filter,  but  tests 
must  be  made  independently  for  each  configuration  because  the  useful 
limit  of  diffusion  will  vary  with  nodal  geometry. 

57 


x  km 


2400 


Figure  16.   48  hour  forecast  for  field  in  figure  15  with  no  diffusion. 
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Figure  17.   Same  as  figure  16  but  with  optimum  diffusion. 
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Figure  18.   Analytic  solution  for  diffusion  tests, 
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Figure  19.   Initial  S  field  for  Robert  filter  tests. 
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C.  EFFECT  OF  ROBERT  FILTER 

For  this  test  the  model  is  configured  nominally  to  maximize  the  accu- 
racy of  the  solution.   The  only  variable  that  is  changed  is  the  Robert 
averaging  coefficient.   Figure  19  is  this  initial  configuration.   The 
averaging  coefficient,  y  starts  at  a  value  of  0.4  and  is  slowly  decreased 
to  0.01  over  11  runs.   The  RMSE  decreases  from  8.7  to  1.74  over  this  range 
with  a  corresponding  increase  in  correlation  and  the  phase  speed  approaches 
unity. 

Since  the  filter  is  desirable  to  control  high  frequency  noise,  a  com- 
promise must  be  reached  that  will  not  unduly  degrade  the  forecast.   For 
the  rest  of  this  investigation  a  coefficient  of  y  =  0.1  is  used  since  the 
associated  RMSE  is  only  2.54.   Figure  20  is  the  48  hour  forecast  field 
with  this  coefficient  and  Figure  21  is  the  corresponding  analytic  solution. 
There  is  little  improvement  in  RMSE  left  to  be  realized  if  any  effective 
filtering  is  desired. 

D.  ACCURACY  VS.  WAVE  NUMBER 

Again  the  model  is  configured  nominally,  except  that  the  number  of 
nodes  is  13  X  24.   Six  runs  are  made  with  the  wave  number  varying  from 
one  to  eight.   One  would  expect  the  forecast  to  deteriorate  as  the  wave 
number  increases.   This  is  due  in  part  to  the  fewer  number  of  nodes  along 
a  latitude  circle  available  to  resolve  each  wave.   But  the  deterioration 
of  RMSE  is  primarily  the  result  of  the  displacement  error  measured  in 
terms  of  the  wavelength.   A  small  displacement  error  with  a  short  wave 
would  produce  a  much  greater  deviation  from  the  true  solution  than  would 
be  produced  by  the  same  displacement  error  in  a  long  wave. 
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Figure  20.   48  hour  forecast  for  field  in  figure  19. 
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Figure  21.   Analytic  solution  for  Robert  filter  tests. 
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The  results  of  this  test  verifies  this  reasoning.   As  the  wave  number 
of  the  initial  field  increases  the  RMSE  worsens,  even  though  the  phase 
speed  error  stays  within  four  percent.   Four  of  these  cases  were  excluded 
from  Section  A  above. 


Phase  Speed 

Wave  Number  Ratio 

1  1.008 

2  1.008 

3  1.02 

4  1.03 

5  1.04 
8  1.03 


E.   RIGHT  TRIANGLES  VS.  EQUILATERAL  TRIANGLES 

Five  pairs  of  runs  are  used  to  test  the  hypothesis  that  elements 
formed  by  equilateral  triangles  yield  forecasts  superior  to  those  based 
upon  elements  formed  by  right  triangles.   In  all  cases  the  flow  is 
uniform  but  the  resolution  is  varied  among  the  cases  from  uniformity 
to  variations  in  either  or  both  directions.   Partial  results  appear  in 
Table  I  as  Cases  1  through  5.   Figures  22,  23  and  24  are  the  initial 
field,  48  hour  forecast  and  analytic  solution,  respectively  for  the 
right  triangles  in  Case  3.   Figures  25,  26  and  27  are  the  corresponding 
figures  for  equilateral  triangles  also  in  Case  3. 

The  differences  between  the  two  arrangements,  while  not  dramatic, 
are  significant  and  consistent.   In  every  case  the  RMSE  is  lower  and 
correlation  higher  for  equilateral  triangles.   The  average  reduction  in 
RMSE  is  about  20  percent.   Not  shown  in  Table  I  but  readily  apparent  in 
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Figure  22.   Initial  S  field  for  Case  3  with  right  triangles, 
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Figure  23.   48  hour  forecast  for  field  in  figure  22, 
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Figure  24.   Analytic  solution  corresponding  to  figure  23, 
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Figure  25.   Initial  S  field  for  Case  3  with  equilateral  triangles. 
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Figure  26.   48  hour  forecast  for  field  in  figure  25 
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Figure  27.   Analytic  solution  corresponding  to  figure  26, 


65 


the  harmonic  analysis,  is  the  generation  of  a  spurious  wave  number  five 
in  the  right  triangle  arrangement  in  four  of  the  five  cases. 

TABLE  I 

RESULTS  OF  RIGHT  VS.  EQUILATERAL  ELEMENT 
AND  DIRECTION  OF  FLOW  TESTS 

RIGHT        EQUILATERAL 


9999 
9996 
9996 
9994 
9992 


CASE 

ri 

r2 

RMSE 

CORR. 

RMSE 

1 

1 

1 

4.25 

.9994 

2.54 

2 

2 

1 

4.54 

.9993 

4.17 

3 

1 

2 

6.40 

.9991 

4.87 

4 

2 

2 

6.96 

.9990 

5.60 

5 

3.55 

1 

6.19 

.9981 

5.75 

6 

1 

4 

7.48 

.9980 

7 

4 

4 

6.60 

.9990 

8 

1 

3.55 

9.00 

.9980 

F.   EFFECT  OF  ADDING  NODES 

Next,  16  runs  are  combined  to  form  nine  test  cases  in  order  to 
determine  the  effect  that  adding  nodes  has  on  the  accuracy  of  the  fore- 
cast.  In  this  section,  the  variables  are  wave  number,  resolution  and 
number  of  nodes.   It  is  found  that  where  the  model  does  very  well  (long 
waves,  uniform  flow,  relatively  small  difference  in  resolution  over  the 
entire  domain) ,  the  addition  of  more  nodes  does  not  increase  the  accuracy, 
especially  if  the  additional  nodes  destroy  the  original  equilateral 
triangularity.   However,  if  the  nodes  are  added  symmetrically  in  both 
the  north-south  and  the  east-west  direction,  so  that  the  effect  is  to 
increase  the  resolution  uniformly  while  preserving  equilateral  triangu- 
larity, then  the  result  is  an  improved  forecast. 
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On  the  other  hand,  when  the  mndel  is  initialized  with  shorter  waves, 
more  nodes  will  improve  the  forecast  even  though  they  may  not  be  added 
symmetrically.   The  major  advantage  of  additional  nodes  is  that  they 
allow  the  resolution  to  be  changed  more  slowly  and  more  smoothly  while 
attaining  a  desired  minimum  resolution.   The  arrangement  with  fewer  nodes 
requires  a  faster  and  more  abrupt  resolution  change  in  order  to  attain 
the  same  minimum  resolution.   The  smoothness  of  the  transition  between 
maximum  and  minimum  resolution  is  apparently  the  critical  factor.   Later 
paragraphs  also  address  this  factor.   The  phase  speed  ratios  for  the 
uniform  grids  are  as  follows: 


Wave 

Phase  Speed 

Nodes 

Number 

Ratio 

RMSE 

13  x  12 

1 

1.001 

2.5 

13  x  24 

1 

1.008 

10.4 

25  x  24  • 

1 

1.003 

7.4 

13  x  12 

2 

1.003 

60.29 

13  x  24 

2 

1.002 

30.7 

25  x  24 

2 

1.003 

25.27 

G.   VARIATION  OF  RESOLUTION  VS.  DIRECTION  OF  FLOW 

The  initial  conditions  for  the  model  are  deliberately  set  to  allow 
only  zonal  flow.   (See  Chapter  V.)   One  of  the  reasons  for  this  is  to 
be  able  to  specifically  test  the  effects  of  variation  of  resolution  along 
the  direction  of  flow.   A  total  of  13  runs  form  six  test  cases.   Many  of 
these  runs  are  the  same  ones  as  used  in  tests  in  preceding  paragraphs. 
Some  of  these  cases  appear  in  Table  I.   In  that  table  the  comparisons 
are  now  vertical.   For  example,  compare  Cases  2  with  Case  3  to  determine 
that  r  =  2  and  r  =  1  is  superior  to  r  =1  and  r  =  2  for  both  right 
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and  equilateral  triangles.   This  means  that  varying  the  resolution  with 
the  flow  is  superior  to  varying  across  it.   Case  2  is  illustrated  in 
Figures  28  and  29  which  are  the  initial  field  and  48  hour  forecast  re- 
spectively.  They  should  be  compared  with  Figures  25  and  26  which,  as 
stated  previously,  are  part  of  Case  3.   The  same  result  holds  for  every 
case  on  the  table.   Additional  tests  have  also  been  run  with  resolution 
ratios  of  six  and  eight.   The  results  are  the  same. 

One  further  test  can  be  conducted  with  the  data  on  Table  I.   Resolu- 
tion ratios  of  r  =  r  =  2  will  yield  the  same  overall  resolution  as  one 
of  the  ratios  equal  to  3.55  and  the  other  unity.   This  is  the  rational 
for  Cases  4,  5  and  8.   These  cases  show  that  better  results  can  be  ob- 
tained by  varying  resolution  in  both  directions  as  opposed  to  only  across 
the  flow  by  an  equivalent  amount.   This  fact  agrees  with  the  results  of 
the  previous  section  in  that  varying  resolution  in  both  directions  tends 
to  maintain  the  equilateral  triangularity  better  than  varying  across  the 
flow  alone. 

H.   SMOOTH  VS.  ABRUPT  VARIATION  OF  RESOLUTION 

For  this  test,  a  version  of  the  model  is  employed  that  allows  the 
resolution  to  be  changed  abruptly.   That  version  is  paired  with  a  more 
standard  version  that  allows  smoother  variation  of  resolution.   Both 
versions  use  right  triangular  elements.   Both  versions  generate  the  same 
maximum  resolution,  but  with  very  different  transitions  to  the  area  of 
coarser  mesh.   Twelve  runs  are  paired  into  six  tests.   Two  different 
wave  numbers  and  three  kinds  of  flow  are  tested.  These  runs  included 
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Figure  28.   Initial  S  field  for  Case  2. 
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Figure  29.   48  hour  forecast  for  Case  2. 
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both  the  normal  and  "flipped"  (See  Section  J  below)  configuration.   The 
results  are  contained  in  Table  II.   Even  though  only  six  tests  are  used, 
the  results  are  dramatic  and  highly  significant.   The  reduction  in  RMSE 
from  the  abrupt  to  the  smooth  cases  averages  62  percent.   Although  not 
shown,  the  reduction  in  bias  averages  70  percent.   Additionally  the 
abrupt  cases  generate  significantly  more  noise  at  all  frequencies.   Cases 
10  and  12  are  the  last  two  cases  excluded  from  Section  A  above.   Figures 
30  through  33  illustrate  the  evolution  of  the  forecast  in  Case  9.   The 
figures  are  the  initial  field  and  12,  24  and  48  hour  forecast  respective- 
ly for  the  abrupt  change.   Compare  these  to  Figures  34  and  35  which  are 
the  initial  field  and  48  hour  forecast  for  the  corresponding  smooth  change, 

Combining  this  result  with  that  of  Section  F  above,  it  is  obvious 
that  the  model  is  highly  sensitive  to  the  rate  and  smoothness  of  the 
transition  from  fine  to  coarse  resolution. 

TABLE  II 
RESULTS  OF  SMOOTH  VS.  ABRUPT  TESTS 


CASE 

9 
10 
11 
12 
11F 
12F 


I.   EFFECT  OF  VARIABLE  WINDS 

In  several  of  the  previous  tests,  reference  has  been  made  to  non- 
uniform flow.   For  this  section  and  the  following  two  sections  a  total 


WAVE 

WIND 

SMOOTH 

ABRUPT 

NUMBER 

RATIO 

RMSE 

CORR. 

RMSE 

CORR 

1 

0.0 

4.54 

.999 

13.6 

.989 

2 

0.0 

28.3 

.956 

89.2 

.470 

1 

0.2 

9.46 

.996 

23.5 

.968 

1 

0.4 

21.4 

.979 

39.9 

.895 

1 

0.2 

4.53 

.999 

14.6 

.987 

1 

0.4 

8.83 

.996 

20.0 

.979 
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Figure  30.   Initial  S  field  with  abrupt  change  of  resolution 
(Case  9)  . 
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Figure  31.   12  hour  forecast  for  field  in  figure  30, 
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Figure  32.   24  hour  forecast  for  field  in  figure  30 
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Figure  33.   48  hour  forecast  for  field  in  figure  30, 
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Figure  34.   Initial  S  field  with  smooth  change  of  resolution 
(Case  9) . 
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Figure  35.   48  hour  forecast  for  field  in  figure  34. 
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of  42  runs  are  made  to  study  the  effects  of  wind  and  the  forcing  term. 
The  42  runs  are  formed  into  14  tests  in  which  the  wind  ratios  are 
specified  as  0.0  for  uniform  flow  and  0.2  and  0.4  for  variable  flow. 
The  tests  include  cases  that  have  both  uniform  and  non-uniform  resolu- 
tion, both  equilateral  and  right  triangular  elements  and  the  13  x  12 
and  13  x  24  nodal  arrangements.    Table  III  contains  the  results  of  12 
of  the  runs  which  comprise  four  tests.   In  12  of  the  total  of  14  tests, 
the  model  performs  better  with  the  uniform  wind.   In  all  tests  once 
the  wind  is  non-uniform  the  forecasts  are  degraded  as  non-uniformity 
increases. 

TABLE  III 
RESULTS  OF  VARIABLE  WIND,  "FLIPPED"  WIND  AND  FORCING  TERM  TESTS 


UNFORCED 

FORCED 

CASE 

RATIO 

RMSE 

SPEED 
. . .NORMAL. . 

CORR. 

RMSE 

SPEED 

CORR. 

13 

0.0 

2.9 

1.004 

.9996 

2.9 

1.004 

.9996 

14 

0.2 

3.4 

1.005 

.9995 

3.3 

1.005 

.9995 

15 

0.4 

8.7 

1.006 
. .FLIPPED. . 

.9970 

8.7 

1.006 

.9970 

16 

0.0 

2.8 

1.005 

.9996 

2.8 

1.007 

.9996 

17 

0.2 

2.7 

1.004 

.9996 

2.6 

1.004 

.9997 

18 

0.4 

3.1 

1.005 

.9995 

2.9 

1.004 

.9960 

J.   "FLIPPED"  WINDS 

Up  until  now,  the  initial  conditions,  in  conjunction  with  the  trans- 
formation of  coordinates  scheme,  require  that  the  area  of  highest  wind 
is  also  the  area  of  highest  resolution,  provided,  of  course,  that  both 
the  wind  and  the  resolution  are  non-uniform.   Now,  downstream  of  the  area 
of  highest  winds  is  an  area  of  convergence  and  downstream  of  the  area 
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of  lowest  wind  is  an  area  of  divergence.   Therefore,  where  the  conver- 
gence runs  into  the  divergence  is  the  area  where  the  gradients  are 
highest.   It  is  a  "bottleneck".   But  that  is  exactly  the  area  where  the 
resolution  is  the  most  coarse.   The  forecast  should  improve  if  the  grid 
or  initial  wind  pattern  could  be  "flipped"  to  concentrate  the  nodes  in 
the  area  of  minimum  winds.   Since  "flipping"  the  grid  requires  fewer 
changes,  that  is  the  choice  here.   As  an  example,  Figures  36  and  37  are 
the  initial  field  and  48  hour  forecast  of  such  an  arrangement. 

To  verify  this  reasoning,  24  runs  are  formed  into  12  tests.   Cases 
13,  14,  and  15  on  Table  III  are  "normal"  whereas  Cases  16,  17  and  18 
are  "flipped".   In  all  tests,  the  "flipped"  arrangement  yields  superior 
forecasts.   This  result  justifies  the  reasoning  that  the  high  resolution 
should  be  concentrated  in  area  of  highest  gradients,  not  the  areas  of 
fastest  flow. 

K.   FORCING  TERM 

The  forcing  term  was  developed  in  Chapter  VIII  as  an  alternate  method 
of  testing  the  forecast  if  the  general  form  of  the  end  product  is  known. 
The  42  runs  form  21  pairs  of  forced/unforced  tests,  some  of  which  are  in 
Table  III.   The  forced  forecasts  are  slightly  better  in  15  tests,  the 
unforced  are  slightly  better  in  five  tests.   One  test  comes  out  a  tie. 
In  general,  the  unforced  forecasts  are  better  if  the  wind  is  uniform  while 
a  non-uniform  wind  favors  the  forced  forecast.   However,  all  differences 
are  very  slight.   Figures  38  and  39  are  the  initial  field  and  48  hour  fore- 
cast respectively  for  a  typical  forced  forecast  with  a  uniform  wind.   They 
should  be  compared  with  Figures  28  and  29  which  form  the  corresponding 

unforced  case. 
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Figure  36.   Initial  S  field  with  "flipped"  winds.   r  =  2,  r  =  1 
and  R  =  0.2.  1       2 
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Figure  37.   48  hour  forecast  for  field  in  figure  36. 
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Figure  38.   Initial  S  field  for  a  typical  forced  case, 
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Figure  39.   48  hour  forecast  for  field  in  figure  38, 
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XI.   CONCLUSIONS 

This  report  applies  the  Finite  Element  Method  to  the  two-dimensional 
advection  equation  with  diffusion.   Although  the  underlying  equation  is 
simple  enough,  FEM  models  have  certain  inherent  advantages,  which  include 
very  accurate  phase  speeds  and  the  ability  to  change  resolution  easily. 
However,  they  also  have  limitations,  such  as  a  relatively  high  computa- 
tional cost,  and  a  requirement  for  fairly  sophisticated  and  efficient 
programming.   The  purpose  of  this  investigation  is  to  establish  some  guide- 
lines and  provide  qualitative  methods  that  may  be  used  in  the  development 
of  future  FEM  models. 

The  conclusions  are  given  below  in  the  same  order  as  they  are  discussed 
in  Chapter  X: 

1.  The  model  yields  reasonable  and  consistent  results.   There  appear  to 
be  no  major  algorithmic  or  coding  problems  as  all  the  results  are 
mathematically  and  meteorologically  plausible. 

2.  Diffusion  is  an  effective  noise  filter;  however  it  does  have  a  maximum 
usable  limit  which  must  be  independently  determined  for  each  applica- 
tion.  Numerical  models  should  contain  at  least  some  provision  for 
diffusion. 

3.  The  Robert  filter  is  also  a  valuable  feature  to  have  built  into  a 
model.   In  general,  the  filter  should  be  no  heavier  than  absolutely 
needed  to  control  high  frequency  noise.   This  model  uses  a  filter  of 
0.1;  higher  values  are  not  recommended. 
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4.  The  model  handles  wave  numbers  1  and  2  well  (24  and  12  nodes  per 
wave,  respectively)  and  wave  number  3  (8  nodes  per  wave)  with  less 
accuracy.   Higher  wave  numbers  maintain  accurate  phase  speeds 
although  RMSE  deteriorates. 

5.  Formulations  with  elements  as  equilateral  triangles  outperform  those 
with  elements  as  right  triangles.   The  difference  in  RMSE  is  about 
20  percent.   The  right  triangle  formulation  generates  spurious  short 
waves  which  may  be  be  the  primary  cause  of  the  relative  inaccuracy. 

6.  Adding  nodes  (and  hence,  elements)  to  the  domain  will  improve  the 
model  if  they  are  added  symmetrically.   More  nodes  allow  smoother 
transition  between  fine  and  coarse  resolution.   This  smoothness  is 
critical. 

7.  An  unexpected  result  is  that  this  model  performs  better  when  the  reso- 
lution is  varied  with  the  flow  rather  than  across  it.   The  reasons 
for  this  are  unclear. 

8.  The  most  important  conclusion  of  this  investigation  is  that  the 
smoother  and  slower  the  change  in  the  resolution,  the  better  the 
forecast.   In  this  model  the  smooth  change  reduces  the  RMSE  by  62 
percent  over  the  abrupt  change. 

9.  The  model  works  best  with  uniform  flow.   In  general,  as  the  flow 
departs  from  uniformity,  the  forecast  continues  to  degrade. 

10.  The  forecast  can  be  improved  by  concentrating  the  high  resolution 
areas  in  the  region  of  the  strongest  gradients. 

11.  Inclusion  of  a  term  that  forces  the  supposed  solution  back  into  the 
model  may  yield  slightly  better  forecasts.   It  forms  a  valid,  alterna- 
tive method  to  test  numerical  models. 
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